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Abstract. There have been some issues in the past in attempts to simulate mag- 
netic fields using the Smoothed Particle Hydrodynamics (SPH) method. SPH is 
well suited to star formation problems because of its Lagrangian nature. We present 
new, stable and conservative methods for magnetohydrodynamics (MHD) in SPH 
and present numerical tests on both waves and shocks in one dimension to show 
that it gives robust and accurate results. 
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1. Introduction 

Smoothed Particle Hydrodynamics (SPH) is a unique method of solving 
the equations of gas dynamics in that it involves no spatial grid. For 
this reason SPH is ideally suited to star formation type problems. 
However there have been some issues in the past with attempts to 
simulate magnetic fields in SPH, namely that when a conservative for- 
mulation of the momentum equation was used an instability developed 
which causes particles to clump together unphysically (Phillips and 
Monaghan, 1985). 

We resolve this instability by adding a short range repulsive force 
to prevent particles from clumping (Monaghan, 2000). In addition we 
formulate the equations of Smoothed Particle Magnetohydrodynam- 
ics (SPMHD) using the continuum equations of Janhunen (2000) and 
Dellar (2001) which are consistent even when the divergence of the 
magnetic field is non zero. Consequently, even though non zero V • B 
may be produced during the simulation, it is treated consistently. 

Shocks are captured within SPMHD by use of artificial dissipation 
terms. The formulation of these terms follows naturally by requiring 
dissipation in the total energy across a shock front and then demanding 
that these terms result in positive definite changes to the entropy. The 
artificial dissipation is effectively turned off away from shocks by use 
of the switch proposed by Morris and Monaghan (1997). 

The resulting equations, when implemented with a simple predictor 
corrector scheme for the timestepping, give good results for a wide range 
of shock tube problems (Price and Monaghan, 2003a; Price and Mon- 
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aghan, 2003b). While we have yet to apply our algorithm to problems 
in two and three dimensions the present results encourage us to believe 
that our SPMHD code will provide a secure basis for astrophysical 
MHD problems. 



2. Numerical tests in one dimension 



The numerical scheme has been tested on a wide variety of one dimen- 
sional problems. We present results here on three standard tests which 
have been used to test many grid-based astrophysical MHD codes (e.g. 
Stone et al., 1992; Dai and Woodward, 1994; Ryu and Jones, 1995; 
Balsara, 1998). For one dimensional MHD problems the magnetic field 
and velocity are allowed to vary in three dimensions. 

The first test involves a fast MHD wave propagating in a periodic 
ID domain and is taken from Dai and Woodward (1998). The wave is 
evolved for 10 periods, corresponding to 10 crossings of the domain. 
The velocity profile in the SPMHD solution is shown in Figure 1, at 
resolutions of 32, 64, 128, 256 and 512 particles (solid points indicate 
the particles). Note that the numerical solutions are in phase with 
the initial conditions (dashed line), demonstrating that the speed at 
which the wave travels is correctly captured by our method. This im- 
provement over previous SPH results (e.g. Morris, 1996) is obtained 
by allowing the smoothing length to vary with local particle density 
and deriving self-consistently (via a Lagrangian variational principle) 
the additional terms which should therefore be included in the SPH 
(Springel and Hernquist, 2002) and correspondingly the SPMHD (Price 
and Monaghan, 2003b) equations. 

The second test (Figure 2) was first described by Brio and Wu 
(1988) and is the MHD analog of the Sod (1978) shock tube problem. 
The problem consists of a discontinuity in pressure, density, transverse 
magnetic field and internal energy initially located at the origin. As 
time develops complex shock structures develop which only occur in 
MHD because of the different wave types. The results shown compare 
well with the solution computed by Balsara (1998) using a grid-based 
MHD code incorporating a Riemann solver, indicated by the solid lines. 

The third example is a test of the code for isothermal MHD and is 
compared with the solution computed by Balsara (1998) in Figure 3. 
The problem illustrates the formation of six discontinuities in the same 
problem. 
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Figure 1. A travelling MHD fast wave propagating in a periodic ID domain. Initial 
conditions are indicated by the dashed line whilst the SPMHD solution is given 
by the solid points after 10 periods (corresponding to 10 crossings of the domain) , 
at resolutions of 32, 64, 128, 256 and 512 particles. The numerical solutions are in 
phase with the initial conditions, demonstrating that the speed at which the wave 
travels is correctly captured by our method. There is a small amount of steepening 
present due to non-linear effects, which is in accordance with the solution presented 
in Dai and Woodward (1998). 



We hope to apply the algorithm to a large simulation of star formation 
similar to that performed by Bate et al. (2003), including the effects 
of magnetic fields. This should provide significant insight into the role 
that magnetic fields play in the star formation process. 
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3. Future work 
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Figure 2. Results of the Brio and Wu (1988) shock tube test. The prob- 
lem consists of an initial discontinuity at the origin, with two phases of gas 
brought into contact at t = 0. The initial conditions to the left of the 
shock are (p, P,v x ,v y , B y ) = [1,1,0,0,1], whilst to the right the conditions are 
(p,P,v x ,v y ,By) = [0.125,0.1,0,0,-1] with B x = 0.75 everywhere and 7 = 2.0. 
These initial conditions are shown by the dashed lines, whilst the solid points, 
corresponding to the SPH particles, indicate the profiles of density, pressure, v x , v y , 
transverse magnetic field B y and thermal energy u at time t = 0.1. The simulation 
uses 800 particles in the one dimensional domain x = [-0.5, 0.5]. The dissipation 
switch has not been used in this case. The solution corresponds well to that computed 
using a grid-based code incorporating a Riemann solver (Balsara 1998 - solid lines). 
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Figure 3. Results of the isothermal MHD shock tube test. This problem illus- 
trates the formation of six discontinuities in isothermal MHD. The setup is 
similar to that in Figure 2 with gas left of the origin initially in the state 
(p,v x ,v y ,v z ,B y ,B z ) = [1.08, 1.2, 0.01, 0.5, 3.6/(4tt) 1/2 , 2/(4tt) 1/2 ] and gas to the 
right in the state (p, v x , v v , v z , B y , B z ) = [1, 0, 0, 0, 4/(4tt) 1/2 , 2/(4tt) 1/2 ] with 
B x — 2/(47r) 1//2 everywhere. These initial conditions are shown by the dashed lines, 
whilst points indicate the positions of the SPH particles at time t = 0.2 which are 
shown to agree well with the numerical solution computed by Balsara (1998) using 
a grid-based code (solid lines). 
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